require(geoR)
## Loading required package: geoR
## --------------------------------------------------------------
##  Analysis of Geostatistical Data
##  For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
##  geoR version 1.9-2 (built on 2022-08-09) is now loaded
## --------------------------------------------------------------
data(s100)
summary(s100)
## Number of data points: 100 
## 
## Coordinates summary
##         Coord.X    Coord.Y
## min 0.005638006 0.01091027
## max 0.983920544 0.99124979
## 
## Distance summary
##         min         max 
## 0.007640962 1.278175109 
## 
## Data summary
##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -1.1676955  0.2729882  1.1045936  0.9307179  1.6101707  2.8678969 
## 
## Other elements in the geodata object
## [1] "cov.model" "nugget"    "cov.pars"  "kappa"     "lambda"
names(s100)
## $coords
## NULL
## 
## $data
## [1] "data"
## 
## $other
## [1] "cov.model" "nugget"    "cov.pars"  "kappa"     "lambda"
var(s100$data)
## [1] 0.7322958
boxplot(s100$data)

## points.geodata()
points(s100)

points(s100, pt.div="equal")

points(s100, cex.min=.6, cex.max=.6)

points(s100, cex.min=.3, cex.max=3)

points(s100, cex.min=1, cex.max=1, pt.div="quart")

points(s100, pt.div="quart")

points(s100, pt.div=4)

points(s100, pt.div="dec")

points(s100, cex.min=1, cex.max=1, pt.div="quint")

points(s100, cex.min=1, cex.max=1, col="gray")

points(s100, cex.min=1, cex.max=1, col=rainbow(100))

points(s100, cex.min=1, cex.max=1, col=terrain.colors(100))

points(s100, pt.div="equal", pch.seq="+")

points(s100, data=sample(s100$data))

points(s100, data=exp(s100$data))

class(s100)
## [1] "geodata"
args(points.geodata)
## function (x, coords = x$coords, data = x$data, data.col = 1, 
##     borders = x$borders, pt.divide = c("data.proportional", "rank.proportional", 
##         "quintiles", "quartiles", "deciles", "equal"), lambda = 1, 
##     trend = "cte", abs.residuals = FALSE, weights.divide = "units.m", 
##     cex.min, cex.max, cex.var, pch.seq, col.seq, add.to.plot = FALSE, 
##     x.leg, y.leg = NULL, dig.leg = 2, round.quantiles = FALSE, 
##     permute = FALSE, ...) 
## NULL
help(points.geodata)
data(parana)
summary(parana)
## Number of data points: 143 
## 
## Coordinates summary
##         east    north
## min 150.1220  70.3600
## max 768.5087 461.9681
## 
## Distance summary
##      min      max 
##   1.0000 619.4925 
## 
## Borders summary
##         east    north
## min 137.9873  46.7695
## max 798.6256 507.9295
## 
## Data summary
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 162.7700 234.1900 269.9200 274.4106 318.2300 413.7000 
## 
## Other elements in the geodata object
## [1] "loci.paper"
borders <- parana$borders
points(parana, border=borders)

points(parana, border=borders, trend="1st")

points(parana, border=borders, trend="1st", pt.div="quart")

## plot.geodata()
plot(s100)

plot(parana)

plot(parana, bor=borders, low=T)

data(ca20)
summary(ca20)
## Number of data points: 178 
## 
## Coordinates summary
##     east north
## min 4957  4829
## max 5961  5720
## 
## Distance summary
##        min        max 
##   43.01163 1138.11774 
## 
## Borders summary
##     east north
## min 4920  4800
## max 5990  5800
## 
## Data summary
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 21.00000 43.00000 50.50000 50.67978 58.00000 78.00000 
## 
## Covariates summary
##     altitude     area   
##  Min.   :3.300   1: 13  
##  1st Qu.:5.200   2: 49  
##  Median :5.650   3:116  
##  Mean   :5.524          
##  3rd Qu.:6.000          
##  Max.   :6.600          
## 
## Other elements in the geodata object
## [1] "reg1" "reg2" "reg3"
names(ca20)
## $coords
## [1] "east"  "north"
## 
## $data
## [1] "data"
## 
## $covariate
## [1] "altitude" "area"    
## 
## $borders
## [1] "borders"
## 
## $other
## [1] "reg1" "reg2" "reg3"
plot(ca20)

points(ca20, cex.min=0.5, cex.max=0.5)
polygon(ca20$reg1)
polygon(ca20$reg2)
polygon(ca20$reg3)

plot(ca20)

plot(ca20,trend=~area)

plot(ca20,trend=~area+altitude)

ksat <- read.geodata("http://www.leg.ufpr.br/geoR/tutorials/datasets/Cruciani.dat", head=T, row.names=1) 
summary(ksat)
## Number of data points: 32 
## 
## Coordinates summary
##        x   y
## min  2.5 0.8
## max 21.0 9.1
## 
## Distance summary
##       min       max 
##  0.509902 18.552897 
## 
## Data summary
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0200000 0.0922500 0.3050000 0.8720625 1.1425000 4.3000000 
## 
## Other elements in the geodata object
## [1] "call"
plot(ksat)

require(MASS)
## Loading required package: MASS
boxcox(ksat)

plot(ksat, lambda=0)

kbor <- read.table("http://www.leg.ufpr.br/geoR/tutorials/datasets/Cruciani.border", head=T, row.names=1) 
kbor
##       x    y
## 1   0.0  0.0
## 2   0.7  1.6
## 3   2.2  3.3
## 4   5.0  5.7
## 5   8.7  7.8
## 6  13.9 10.4
## 7  19.8 13.0
## 8  21.0 10.5
## 9  22.4  5.3
## 10 22.5  0.0
plot(ksat, lambda=0, bor=kbor)

points(ksat, lambda=0, bor=kbor)

## Variog
s100.v <- variog(s100, max.d=1)
## variog: computing omnidirectional variogram
plot(s100.v)

parana.v <- variog(parana, max.dist=400, trend="1st")
## variog: computing omnidirectional variogram
plot(parana.v)

ca20.v <- variog(ca20, trend=~area, max.d=800)
## variog: computing omnidirectional variogram
plot(ca20.v)

ksat.v <- variog(ksat, max.d=10, lam=0)
## variog: computing omnidirectional variogram
plot(ksat.v)